Fit example models

Author

Max Lindmark

Published

2024-03-28

# Load libraries
library(tidyverse)
library(tidylog)
library(sdmTMB)
library(patchwork)
library(viridis)
library(RColorBrewer)
library(ggsidekick); theme_set(theme_sleek())

home <- here::here()
# source_url("https://raw.githubusercontent.com/maxlindmark/cod-interactions/main/R/functions/map-plot.R")
source(paste0(home, "/R/functions/map-plot.R"))

Read data and prediction grid, scales variables

# Read data
d <- readr::read_csv(paste0(home, "/data/clean/larval_size.csv")) %>% 
  drop_na(temp) %>% 
  mutate(temp_sc = (temp - mean(temp, na.rm = TRUE)) / sd(temp, na.rm = TRUE),
         yday_ct = yday - mean(yday),
         year_f = as.factor(year))
Rows: 39674 Columns: 14
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (3): haul_id, species, period
dbl (11): length_mm, year, day, month, lat, lon, temp_obs, yday, X, Y, temp

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
drop_na: removed 96 rows (<1%), 39,578 rows remaining

mutate: new variable 'temp_sc' (double) with 1,488 unique values and 0% NA

        new variable 'yday_ct' (double) with 46 unique values and 0% NA

        new variable 'year_f' (factor) with 28 unique values and 0% NA
nrow(d)
[1] 39578
pred_grid <- readr::read_csv(paste0(home, "/data/clean/pred_grid.csv")) %>% 
  drop_na(temp) %>% 
  mutate(temp_sc = (temp - mean(d$temp, na.rm = TRUE)) / sd(d$temp, na.rm = TRUE),
         year_f = as.factor(year),
         yday_ct = 0)
Rows: 144971 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (7): X, Y, year, lon, lat, depth, temp

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
drop_na: removed 1,972 rows (1%), 142,999 rows remaining

mutate: new variable 'temp_sc' (double) with 125,655 unique values and 0% NA

        new variable 'year_f' (factor) with 29 unique values and 0% NA

        new variable 'yday_ct' (double) with one unique value and 0% NA

Fit example model

pred_grid_list <- list()
coef_list <- list()
res_list <- list()
cond_list <- list()

for(i in unique(d$species)) {  
    
    dd <- d %>% filter(species == i)
    
    mesh <- make_mesh(dd,
                      xy_cols = c("X", "Y"),
                      cutoff = 5)
    
    print(
    ggplot() +
      inlabru::gg(mesh$mesh) +
      coord_fixed() +
      geom_point(aes(X, Y), data = dd, alpha = 0.2, size = 0.5) +
      annotate("text", -Inf, Inf, label = paste("n knots = ", mesh$mesh$n), hjust = -0.1, vjust = 2) + 
      labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(i, "_", " ")))
    )
 
    # TODO: Time-varying intercept (rw or ar1) to fill in gaps with ar1?
    # TODO: Anisotrophy might be warranted here
    # TODO: Spatiotemporal?
    # TODO: Spatially varying NAO?
    # TODO: R2 for temperature (write a functions)
    
    if(unique(!dd$species %in% c("Maurolicus muelleri",
                                 "Anguilla anguilla",
                                 "Agonus cataphractus",
                                 "Syngnathus rostellatus",
                                 "Enchelyopus cimbrius"))) {
      
    m <- sdmTMB(length_mm ~ 0 + year_f + temp_sc + yday_ct,
                data = dd,
                mesh = mesh,
                family = student(df = 8),
                spatiotemporal = "off",
                spatial = "on",
                time = "year")  
    
    # Residuals
    samps <- sdmTMBextra::predict_mle_mcmc(m, mcmc_iter = 201, mcmc_warmup = 200)
    dd$mcmc_res <- as.vector(residuals(m, type = "mle-mcmc", mcmc_samples = samps))
      
    } else {
      
    m <- sdmTMB(length_mm ~ 0 + year_f + temp_sc + yday_ct,
                data = dd,
                family = student(df = 8),
                spatiotemporal = "off",
                spatial = "off",
                time = "year")  
      
    # Residuals
    # FIXME: why doesn't MCMC residuals work without spatial models?
    dd$mcmc_res <- residuals(m)
    
    }
    
    print(i)
    m
    sanity(m)
    
  
    res_list[[i]] <- dd %>% mutate(species = i)
    
    # Gridded predictions
    pred <- predict(m, newdata = pred_grid %>% filter(year %in% unique(dd$year)))
    
    pred_grid_list[[i]] <- pred %>% mutate(species = i)

    # Coefficients
    coef_list[[i]] <- tidy(m, conf.int = TRUE) %>%
      filter(term %in% c("temp_sc", "yday_ct")) %>% 
      mutate(species = i)
    
    # Conditional predictions
    nd <- tibble(year = unique(dd$year)) %>% 
      mutate(year_f = as.factor(year), 
             yday_ct = 0, 
             temp_sc = 0)
    
    pred_c <- predict(m, newdata = nd, se_fit = TRUE,
                      re_form = NA, re_form_iid = NA)

    cond_list[[i]] <- pred_c %>% mutate(species = i)
    
}
filter: removed 38,358 rows (97%), 1,220 rows remaining
Loading required namespace: INLA

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.001667 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 16.67 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 5.052 seconds (Warm-up)
Chain 1:                0.017 seconds (Sampling)
Chain 1:                5.069 seconds (Total)
Chain 1: 
[1] "Ammodytidae"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
mutate: no changes
filter: removed 4,931 rows (3%), 138,068 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 28 rows (93%), 2 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'year_f' (factor) with 28 unique values and 0% NA
        new variable 'yday_ct' (double) with one unique value and 0% NA
        new variable 'temp_sc' (double) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 29,496 rows (75%), 10,082 rows remaining


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.003169 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 31.69 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 13.447 seconds (Warm-up)
Chain 1:                0.032 seconds (Sampling)
Chain 1:                13.479 seconds (Total)
Chain 1: 
[1] "Clupea harengus"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
mutate: no changes
filter: removed 4,931 rows (3%), 138,068 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 28 rows (93%), 2 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'year_f' (factor) with 28 unique values and 0% NA
        new variable 'yday_ct' (double) with one unique value and 0% NA
        new variable 'temp_sc' (double) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 33,778 rows (85%), 5,800 rows remaining


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.002786 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 27.86 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 10.101 seconds (Warm-up)
Chain 1:                0.031 seconds (Sampling)
Chain 1:                10.132 seconds (Total)
Chain 1: 
[1] "Aphia minuta"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
mutate: no changes
filter: removed 4,931 rows (3%), 138,068 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 28 rows (93%), 2 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'year_f' (factor) with 28 unique values and 0% NA
        new variable 'yday_ct' (double) with one unique value and 0% NA
        new variable 'temp_sc' (double) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 38,970 rows (98%), 608 rows remaining


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.001133 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 11.33 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 17.636 seconds (Warm-up)
Chain 1:                0.043 seconds (Sampling)
Chain 1:                17.679 seconds (Total)
Chain 1: 
[1] "Microstomus kitt"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
mutate: no changes
filter: removed 19,724 rows (14%), 123,275 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 25 rows (93%), 2 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'year_f' (factor) with 25 unique values and 0% NA
        new variable 'yday_ct' (double) with one unique value and 0% NA
        new variable 'temp_sc' (double) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 28,532 rows (72%), 11,046 rows remaining


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.003882 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 38.82 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 13.538 seconds (Warm-up)
Chain 1:                0.045 seconds (Sampling)
Chain 1:                13.583 seconds (Total)
Chain 1: 
[1] "Crystallogobius linearis"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
mutate: no changes
filter: removed 4,931 rows (3%), 138,068 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 28 rows (93%), 2 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'year_f' (factor) with 28 unique values and 0% NA
        new variable 'yday_ct' (double) with one unique value and 0% NA
        new variable 'temp_sc' (double) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 39,523 rows (>99%), 55 rows remaining

[1] "Maurolicus muelleri"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
mutate: no changes
filter: removed 93,689 rows (66%), 49,310 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 10 rows (83%), 2 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'year_f' (factor) with 10 unique values and 0% NA
        new variable 'yday_ct' (double) with one unique value and 0% NA
        new variable 'temp_sc' (double) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 38,503 rows (97%), 1,075 rows remaining


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.001179 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 11.79 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 2.734 seconds (Warm-up)
Chain 1:                0.012 seconds (Sampling)
Chain 1:                2.746 seconds (Total)
Chain 1: 
[1] "Chirolophis  ascanii"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
mutate: no changes
filter: removed 4,931 rows (3%), 138,068 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 28 rows (93%), 2 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'year_f' (factor) with 28 unique values and 0% NA
        new variable 'yday_ct' (double) with one unique value and 0% NA
        new variable 'temp_sc' (double) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 37,648 rows (95%), 1,930 rows remaining


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.001625 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 16.25 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 10.658 seconds (Warm-up)
Chain 1:                0.036 seconds (Sampling)
Chain 1:                10.694 seconds (Total)
Chain 1: 
[1] "Pomatoschistus sp"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
mutate: no changes
filter: removed 9,862 rows (7%), 133,137 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 27 rows (93%), 2 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'year_f' (factor) with 27 unique values and 0% NA
        new variable 'yday_ct' (double) with one unique value and 0% NA
        new variable 'temp_sc' (double) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 39,311 rows (99%), 267 rows remaining

[1] "Anguilla anguilla"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
mutate: no changes
filter: removed 54,241 rows (38%), 88,758 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 18 rows (90%), 2 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'year_f' (factor) with 18 unique values and 0% NA
        new variable 'yday_ct' (double) with one unique value and 0% NA
        new variable 'temp_sc' (double) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 35,409 rows (89%), 4,169 rows remaining


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.001747 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 17.47 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 5.218 seconds (Warm-up)
Chain 1:                0.018 seconds (Sampling)
Chain 1:                5.236 seconds (Total)
Chain 1: 
[1] "Pholis gunnellus"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
mutate: no changes
filter: removed 4,931 rows (3%), 138,068 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 28 rows (93%), 2 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'year_f' (factor) with 28 unique values and 0% NA
        new variable 'yday_ct' (double) with one unique value and 0% NA
        new variable 'temp_sc' (double) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 39,199 rows (99%), 379 rows remaining

[1] "Agonus cataphractus"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
mutate: no changes
filter: removed 34,517 rows (24%), 108,482 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 22 rows (92%), 2 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'year_f' (factor) with 22 unique values and 0% NA
        new variable 'yday_ct' (double) with one unique value and 0% NA
        new variable 'temp_sc' (double) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 38,466 rows (97%), 1,112 rows remaining

[1] "Syngnathus rostellatus"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
mutate: no changes
filter: removed 14,793 rows (10%), 128,206 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 26 rows (93%), 2 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'year_f' (factor) with 26 unique values and 0% NA
        new variable 'yday_ct' (double) with one unique value and 0% NA
        new variable 'temp_sc' (double) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 38,650 rows (98%), 928 rows remaining


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000756 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 7.56 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 5.344 seconds (Warm-up)
Chain 1:                0.013 seconds (Sampling)
Chain 1:                5.357 seconds (Total)
Chain 1: 
[1] "Sprattus sprattus"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
mutate: no changes
filter: removed 59,172 rows (41%), 83,827 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 17 rows (89%), 2 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'year_f' (factor) with 17 unique values and 0% NA
        new variable 'yday_ct' (double) with one unique value and 0% NA
        new variable 'temp_sc' (double) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 39,502 rows (>99%), 76 rows remaining


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000205 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.05 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 2.494 seconds (Warm-up)
Chain 1:                0.006 seconds (Sampling)
Chain 1:                2.5 seconds (Total)
Chain 1: 
[1] "Argentina spp"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
mutate: no changes
filter: removed 98,620 rows (69%), 44,379 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 9 rows (82%), 2 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'year_f' (factor) with 9 unique values and 0% NA
        new variable 'yday_ct' (double) with one unique value and 0% NA
        new variable 'temp_sc' (double) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 39,070 rows (99%), 508 rows remaining


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000661 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 6.61 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 2.615 seconds (Warm-up)
Chain 1:                0.012 seconds (Sampling)
Chain 1:                2.627 seconds (Total)
Chain 1: 
[1] "Myoxocephalus scorpius"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
mutate: no changes
filter: removed 29,586 rows (21%), 113,413 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 23 rows (92%), 2 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'year_f' (factor) with 23 unique values and 0% NA
        new variable 'yday_ct' (double) with one unique value and 0% NA
        new variable 'temp_sc' (double) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 39,510 rows (>99%), 68 rows remaining

[1] "Enchelyopus cimbrius"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
mutate: no changes
filter: removed 98,620 rows (69%), 44,379 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 9 rows (82%), 2 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'year_f' (factor) with 9 unique values and 0% NA
        new variable 'yday_ct' (double) with one unique value and 0% NA
        new variable 'temp_sc' (double) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 39,323 rows (99%), 255 rows remaining


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000306 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 3.06 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 2.958 seconds (Warm-up)
Chain 1:                0.005 seconds (Sampling)
Chain 1:                2.963 seconds (Total)
Chain 1: 
[1] "Limanda limanda"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
mutate: no changes
filter: removed 93,689 rows (66%), 49,310 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
filter: removed 10 rows (83%), 2 rows remaining
mutate: new variable 'species' (character) with one unique value and 0% NA
mutate: new variable 'year_f' (factor) with 10 unique values and 0% NA
        new variable 'yday_ct' (double) with one unique value and 0% NA
        new variable 'temp_sc' (double) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA

# Save predictions and sims as data frames
res_df <- dplyr::bind_rows(res_list)
pred_grid_df <- dplyr::bind_rows(pred_grid_list)
coef_df <- dplyr::bind_rows(coef_list)
cond_df <- dplyr::bind_rows(cond_list)

Save

write_csv(res_df, paste0(home, "/output/res_df.csv"))
write_csv(pred_grid_df, paste0(home, "/output/pred_grid_df.csv"))
write_csv(coef_df, paste0(home, "/output/coef_df.csv"))
write_csv(cond_df, paste0(home, "/output/cond_df.csv"))

Plot output

res_df <- read_csv(paste0(home, "/output/res_df.csv"))
pred_grid_df <- read_csv(paste0(home, "/output/pred_grid_df.csv"))
coef_df <- read_csv(paste0(home, "/output/coef_df.csv"))
cond_df <- read_csv(paste0(home, "/output/cond_df.csv"))

Residuals

res_df %>% 
  ggplot(aes(sample = mcmc_res)) +
  stat_qq(size = 0.75, shape = 21, fill = NA) +
  facet_wrap(~species) +
  stat_qq_line() +
  labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
  theme(aspect.ratio = 1)

ggsave(paste0(home, "/figures/residuals.pdf"), width = 17, height = 17, units = "cm")

Coefficients

pal <- brewer.pal(n = 6, name = "Paired")[c(2, 6)]

coef_df %>% 
  filter(term == "temp_sc") %>% 
  mutate(sign = ifelse(estimate > 0, "pos", "neg"),
         sig = ifelse(estimate > 0 & conf.low > 0, "sig", "not sig"),
         sig = ifelse(estimate < 0 & conf.high < 0, "sig", sig)) %>% 
  ggplot(aes(estimate, reorder(species, estimate), fill = sig, color = sign, shape = sig)) + 
  geom_point(fill = NA) + 
  scale_shape_manual(values = c(21, 19)) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), 
                width = 0, alpha = 0.3) + 
  geom_vline(xintercept = 0, alpha = 0.3, linetype = 2) +
  theme(axis.text.y = element_text(face = "italic")) + 
  labs(y = "Species", x = "Effect of temperature") + 
  scale_color_manual(values = pal) + 
  scale_fill_manual(values = pal) +
  guides(color = "none", shape = "none") +
  theme(legend.position = "bottom")

ggsave(paste0(home, "/figures/temp_coef.pdf"), width = 17, height = 8, units = "cm")

Random fields

pred_grid_df_sub <- pred_grid_df %>% 
  group_by(species) %>% 
  filter(year == max(year))

plot_map_fc +
  geom_raster(data = pred_grid_df_sub %>% drop_na(omega_s),
              aes(X*1000, Y*1000, fill = omega_s)) +
  scale_fill_gradient2() + 
  theme_facet_map(base_size = 8) +
  facet_wrap(~species, ncol = 4)

ggsave(paste0(home, "/figures/omega.pdf"), width = 17, height = 15, units = "cm")

Spatial predictions

For all years

# Z- score sizes by species to plot on the same plot

plot_map_fc +
  geom_raster(data = pred_grid_df %>% 
                drop_na(omega_s) %>% # Drop the species without spatial effects
                mutate(est_z = (est - mean(est)) / sd(est), .by = "species"), 
              aes(X*1000, Y*1000, fill = est_z)) +
  scale_fill_viridis() + 
  theme_facet_map(base_size = 8) +
  facet_wrap(~species, ncol = 4)

ggsave(paste0(home, "/figures/spatial_pred.pdf"), width = 17, height = 15, units = "cm")